Modified spin-wave theory with ordering vector optimization I: frustrated 
bosons on the spatially anisotropic triangular lattice 

Philipp Hauke, 1,2, [^Tommaso Roscilde, 3 Valentin Murg, 2 J. Ignacio Cirac, 2 and Roman Schmied 2 

1 ICFO - Institut de Ciencies Fotdniques, 
Av. Canal OUmpic s/n, E-08860 Castelldefels (Barcelona), Spain 
2 Max-Planck-Institut fur Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany 
3 Laboratoire de Physique, Ecole Normale Superieure de Lyon, 46 Allee d'ltalie, F-69007 Lyon, France 

(Dated: December 29, 2009) 

We investigate a system of frustrated hardcore bosons, modeled by an XY antiferromagnet 
on the spatially anisotropic triangular lattice, using Takahashi's modified spin-wave (MSW) 
theory. In particular we implement ordering vector optimization on the ordered reference 
state of MSW theory, which leads to significant improvement of the theory and accounts for 
quantum corrections to the classically ordered state. The MSW results at zero temperature 
compare favorably to exact diagonalization (ED) and projected entangled-pair state (PEPS) 
calculations. The resulting zero-temperature phase diagram includes a ID quasi-ordered 
phase, a 2D Neel ordered phase, and a 2D spiraling ordered phase. Strong indications com- 
ing from the ED and PEPS calculations, as well as from the breakdown of MSW theory, 
suggest that the various ordered or quasi-ordered phases are separated by spin-liquid phases 
with short-range correlations, in analogy to what has been predicted for the Heisenberg 
model on the same lattice. Within MSW theory we also explore the finite-temperature 
phase diagram. In agreement with Bcrczinskii-Kosterlitz- Thouless (BKT) theory, we find 
that zero-temperature long-range-ordered phases turn into quasi-ordered phases (up to a 
BKT transition temperature), while zero-temperature quasi-ordered phases become short- 
range correlated at finite temperature. These results show that, despite its simplicity, mod- 
ified spin-wave theory is very well suited for describing ordered and quasi-ordered phases 
of frustrated XY spins (or, equivalently, of frustrated lattice bosons) both at zero and fi- 
nite temperatures. While MSW theory, just as other theoretical methods, cannot describe 
spin-liquid phases, its breakdown provides a fast and reliable method for singling out Hamil- 
tonians which may feature these intriguing quantum phases. We thus suggest a tool for 
guiding our search for interesting systems whose properties are necessarily studied with a 
physical quantum simulator instead of theoretical methods. 
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I. INTRODUCTION 

Lattice models of strongly interacting bosons have recently been implemented experimentally 
thanks to impressive developments with trapped ultra-cold atoms in optical lattice potentials [U 
[2]. A particularly appealing perspective in this arena is the study of strongly correlated lattice 
boson models with frustration, arising for instance from the coupling of bosons to a (artificial) 
magnetic field HI or from a periodical shaking of the optical lattice [6 ]. Frustration in 
the intersite hopping amplitudes is formally equivalent to a description of the system in a rotating 
reference frame, which implies that the system is subject to the spontaneous appearance of vortices. 
Such vortices can form ordered arrays (vortex crystals) coexisting with Bose condensation, which 
consequently takes place in a macroscopic wavefunction sustaining persisting circulating currents 
(see Ref . [7] and references therein) ; or they can even disrupt condensation completely, and lead to a 
disordered insulating state [8 j . Such disordered states are notoriously difficult to study theoretically. 

In the particular limit of a very strong interparticle repulsion, frustrated bosonic models can 
be exactly mapped onto S = 1/2 frustrated XY antiferromagnets [9]. In two dimensions these 
models are known to exhibit ground states with spiral order, representing the magnetic counterpart 
to the aforementioned Bose-condensed states with vortex arrays. In special circumstances the 
interplay between quantum fluctuations and frustration may lead to disordered spin-liquid states, 
which are in one-to-one correspondence with bosonic insulating phases. XY antiferromagnets can 
also be regarded as the limiting case of antiferromagnetic Hamiltonians with planar anisotropy in 
the couplings, relevant to the description of frustrated antiferromagnetic materials, and they can 
describe the physics of Cooper pairs in arrays of ultra-small Josephson junctions with magnetic 
frustration |1U] . More recently we have proposed that frustrated XY antiferromagnets can be 
experimentally implemented by loading planarly trapped ions into an optical lattice jllj . 

From a theoretical point of view, bosonic frustration in the presence of strong interparticle 
interactions on a lattice represents a very hard problem in dimensions d > 1, due to the lack of 
controlled perturbative expansions in the strongly correlated regime; to the breakdown of semi- 
classical methods in the presence of strong quantum fluctuations enhanced by frustration; and to 
the appearance of a sign problem in quantum Monte Carlo simulations. Hence the implementation 
of bosonic frustration in optical lattice experiments would represent a fundamental instance of a 
useful quantum simulation, possibly outperforming any classical computation (see e.g. Ref. |12j 
and references therein). Indeed, as mentioned above, solving frustrated bosonic models amounts 
to solving a large class of frustrated antiferromagnets. Fundamental steps have very recently been 
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taken experimentally towards the implementation of artificial magnetic fields in cold atom experi- 
ments via Raman schemes [13]. Hence exciting progress in this field is expected in the near future. 
However, the difficulty of finding accurate theoretical descriptions of disordered quantum lattice 
models makes it hard to tell a priori which systems will present such interesting phases in an 
experiment. The most attractive aspect of quantum simulators is their potential ability to emulate 
model Hamiltonians whose phase diagram cannot be accurately predicted with current theoreti- 
cal approaches. Therefore it would be highly desirable to dispose of a fast tool that can outline 
quantum-mechanical phase diagrams, point out disordered phases, and thus classify model Hamil- 
tonians according to their potential interest for experimental quantum simulation. We propose 
that the methods presented here will serve this very purpose. 

Planar systems of bosons in optical lattices, coupled to an artificial magnetic field, can be 
described by the Bose-Hubbard Hamiltonian 

^bh = £ \ (bib, + h.c.) + | J2 n ^ ~ !) (!) 

(hi) * 

where bi, b\ are bosonic operators, rtj = b\b{, and represents pairs of nearest neighbor sites. 
The hopping amplitude ty = —tij exp (iAij) , where Uj > 0, is spatially modulated by the line 
integral of the vector potential along the bond, Aij = JJ* 1 A(r) ■ dl, as well as by possible 
spatial anisotropies in the optical lattice (contained in the bare hopping amplitudes tij > 0). In 
the following we will focus on the limit of infinite repulsion U — ► oo and half filling (m) = 1/2, 
under which the Bose-Hubbard model maps onto the S = 1/2 XY Hamiltonian |llj : 

V, Y.'u(^"Sf ■ Sf'Sf) (2) 

(id) 

where S-* is the a th component of the S = 1/2 spin operator acting on site i. 

In this work we focus on a triangular lattice with antiferromagnetic nearest-neighbor interac- 
tions, which can be seen as a ferromagnetic lattice (tij > 0) with half a magnetic flux quantum 
threaded through each lattice plaquette. This magnetic flux can, for instance, be interpreted as 
flipping the signs of all hopping amplitudes of the bonds along the horizontal direction of the lattice 
(Fig. 1). In this gauge the hopping is transformed to = tij for r{ — rj = ±<zti, where a is the 
lattice spacing and t\ = (1,0), while tij = —tij for the other (diagonal) nearest- neighbor bond 
directions. A canonical transformation rotating the spins in odd rows by an angle tt around the z 
axis flips the signs of all diagonal bonds, leading to a model of a spatially anisotropic triangular 
XY antiferromagnet with all tij > 0. Such a model exhibits strong frustration on the triangular 
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Figure 1: Spatially anisotropic triangular lattice. The lines denote interactions and the spins are located at 
the vertices. 

lattice. In the following we specialize to the case in which the hopping amplitudes ty take two 
values, ti and ti depending on the orientation of the bond [parallel to t\ or along the diagonals 
±ar2 and ±a (T2 — Ti), with T2 = (l/2, \/3/2) - see Fig.jlj. This particular case is relevant for the 
physics of a number of systems: it can describe neutral atoms trapped in triangular lattices formed 
with three lasers intersecting at 120° angles [14 , and one of which has a different intensity than 
the other two; another implementation of this model are neutral atoms in an isotropic triangular 
lattice in the absence of artificial magnetic fields but with elliptical time-dependent forcing of the 
lattice [6] ; furthermore, it is relevant for triangular Wigner crystals of ions trapped in the minima 
of an optical lattice . 

The aim of this work is the determination of the ground-state and finite-temperature phase 
diagrams of the S = 1/2 XY antiferromagnet (AF) - or, alternatively, of frustrated half-filled 
hardcore bosons - on a spatially anisotropic triangular lattice (SATL) by means of spin-wave 
theory. In particular, we show that Takahashi's modified spin-wave (MSW) theory [15^ gives an 
adequate description of the main features of the ground-state and low-temperature phase diagram, 
while keeping the computational effort at a minimum. Spin-wave methods generally account for 
weak quantum fluctuations around the ordered state corresponding to the ground state in the 
classical limit S — > 00. We show how MSW theory can be extended to arbitrary reference states, 
whose ordering vector may be shifted with respect to the ground state in the classical limit |16j . 
Furthermore, we demonstrate how this procedure allows for a convenient calculation of the spin 
stiffness. In the specific case of the S = 1/2 XYAF on a SATL, the spin stiffness proves to be 
an effective tool supporting our search for spin-liquid phases. Supplementary results derived by 
projected entangled-pair states (PEPS) [17 and exact diagonalizations (ED) using the Lanczos 
method allow us to validate the zero-temperature results of the MSW method with ordering vector 
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optimization. We moreover provide the finite-temperature phase diagram of the S = 1/2 XYAS 
on a SATL, and find a region where the breakdown of MSW theory indicates a disordered state. 
Such a disordered phase is an ideal candidate for performing a useful quantum simulation because 
theoretical tools for studying it adequately are currently lacking. In an upcoming publication we 
will extend this method to the ground state phase diagrams of the Heisenberg SATL and the J1J2J3 
model HBJ. 

For coherence with the theoretical technique used - spin-wave theory - and for a better com- 
parison with existing results in the literature, the results of this paper will be generally expressed 
in the language of spin physics, but guidance will be provided on how to translate the magnetic 
observables into bosonic observables. 

This paper is organized as follows. In section [IT] we introduce Takahashi's MSW formalism, 
supplemented with ordering-vector optimization, and provide a general method to calculate the 
spin stiffness. The rest of the work is dedicated to the investigation of the phase diagram of S = 1/2 



XYAF on a SATL. Section III presents the ground state phase diagram of this model. In section IV 



we extend the phase diagram to finite temperatures and calculate Berezinskii-Kosterlitz-Thouless 
transition lines. Finally, in section IV] we present our conclusions. 



II. MODIFIED SPIN- WAVE FORMALISM 



In the past MSW theory has been found to give a satisfactory qualitative account of the low- 
temperature properties of spin systems, even frustrated or disordered ones. In this section we 
review its formalism, mainly following Xu and Ting [16] in the first steps, but considering XY 
interactions rather than Heisenberg interactions (for the latter see also Ref. [IS). This requires 
only minor modifications of the formulas, and it is expected that the validity of the spin-wave 
approach is even better justified for XY interactions since in this case the influence of quantum 
fluctuations is reduced by the anisotropy in spin coupling. 

Our aim is to determine the phase diagram of the Hamiltonian of Eq. (|2]). A fundamental 
assumption of spin-wave theory as applied to the XY model is that the ground state shows long- 
range order (LRO) with the spins classically lying in the xy-plane; for a translationally invariant 
system, like the one under investigation, the ordered ground state is characterized by a well defined 
ordering vector Q. Under this assumption it is convenient to rotate the local reference system of 
each spin from {x, y, z) to (77, £) so that the ground state in the local reference frame has all spins 
aligned in the same direction. This amounts to the following transformation: 
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Sf = -sm(Q- ri )SP + cos (Q-ri)S> , (3a) 
S? = cos(Q-r i )5/ ? +sin(Q-r l )S/, (3b) 
Si = -Sf . (3c) 

Then & , which will be the quantization axis, lies parallel to the classical spin Si = 
(cos (Q ■ Ti) ,sin (Q • r%) , 0). The component lies perpendicular to it in the xy plane, and 

is perpendicular to the xy plane. Unlike in ordinary spin- wave theories we do not make any 
assumption on the ordering vector Q. In particular, it may well differ from the one exhibited by 
the system in the classical limit (Q cl ). 

Spin waves can be described by applying the Dyson-Maleev (DM) transformation |19| [2"U] , which 
maps the physical spins to interacting bosons, 

Sr - -^fa-ala^Oi, (4a) 
S+ -» V2Sa\, (4b) 
5/ -> -S + alai, (4c) 

where = s£ ± iS? '■ The DM transformation is an exact mapping between spins and bosons as 
long as projectors are retained which keep the system in the physical subspace, i.e. the subspace 
where at each site only 2S DM bosons are present at most. It can be shown that these projectors 
have the form V = 1 + 0[n/(2S)] 3 where n is the DM boson density |21j . Hence, under the 
assumption of diluteness of the DM boson gas, n/(2S) < 1 (in fact (n) = S, see below), we can 
safely neglect the V projectors. 

If the spin Hamiltonian under investigation is obtained as the hardcore limit of the Bose- 
Hubbard Hamiltonian Eq. ([T]), it is important to distinguish the DM bosons from the physical 
6-bosons from which the effective spin Hamiltonian originated. Indeed the DM bosons at a site i 
quantify the deviation of the i th spin from the local direction in the xy plane set by the ordered 
structure with ordering vector Q. On the other hand, the physical bosons correspond to the spin 
deviations with respect to full alignement of the spin along the z axis. The particle- hole symmetry 
of the bosonic Hamiltonian, Eq. ([!]), in the hardcore limit leads to half filling of the physical bosons, 
which accidentally coincides with the average filling imposed by the Takahashi's constraint on the 
DM bosons for 5 = 1/2 (see next subsection). Yet all other properties are in general quite different. 
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A. Derivation of a mean-field Hamiltonian and Takahashi's constraint 



Applying Eqs. and @ to Eq. ^ one arrives at the bosonic Hamiltonian 



H = S 2 ^2tij 1 - — ^2ajaj + 2aja,- - ajaj - a^a] + a\a^ + 



+ 



(2S) 2 \ a \ aia ) a i ~ a i a ) a i a i ~ a i a i a i a j + aia^ajCLj + ajajdja.,-J + O 



,t„t. 



/ n \ 3 



cos (Q • Tit) . (5) 



~j~J — i— *— *— j i 1 j 

Here we have dropped the terms with six boson operators, which are of order 0[n/(25) 3 ] and are 
negligible for n/{2S) < 1. Moreover the truncation of the Hamiltonian to this order is consis- 
tent with neglecting the effect of the projectors on the physical subspace which also amounts to 
discarding terms of order 0[n/(2S) 3 ]. 

MSW theory relies on the minimization of the free energy. To this end we need the expectation 
value of the Hamiltonian Eq. Under the assumption that the ground state is a Gaussian state 
we make use of Wick's theorem |22j to decouple the boson-boson interaction terms, i.e. the terms 
consisting of four boson operators. The expectation value E = (TC) can then be written as 



1 \ - 



S+ 2 -F(0)+F( rij ) 



+ 



S+ 2 -F(0) + G( riJ ) 



COS (Q ■ Ti 



Here we have defined the correlators 



(4 a j) = F ( r ij) ~ 

(a iaj ) = (ojaj) = G(r«)- 



(6) 

(7a) 
(7b) 



"01 ~ V**"^/ ~ w V* ijj 

These correlators can be rewritten in terms of independent particles by first Fourier transforming 
a/, = — ^ Yli a i eT lk ' Vi , where N is the number of sites, and then applying a Bogoliubov transfor- 



mation 



Oik 



cosh 8k dk — sinh Ok a_ k 



n 



sinh Ok afc + cosh Ok a 



Requiring that the Bogoliubov particles be non-interacting imposes (afcOfc' 
condition also removes the anti-Hermitian part of the Hamiltonian. 
The correlators are now 



(44 



(8a) 
(8b) 

0. This 



(9a) 
(9b) 
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with n k = (a k a k ) = 1/ (exp(u;fc/T) — 1) being the occupation number of Bogoliubov mode k 
at temperature T (with the Boltzmann constant k-g set to unity). The dispersion relation uo k is 
determined self-consistently in the next section. 

So far we have essentially formulated a standard Hartree-Fock theory for the gas of interacting 
DM bosons. A very important modification to this theory, due to Takahashi [T5], is the introduction 
of the constraint of zero magnetization at each site, 

<S/> = -S + {a\ ai ) = -S - \ + F (0) = 0. (10) 

The implementation of this constraint amounts to effectively reducing the Hilbert space dimension 
available to the DM bosons by fixing their average density to S. For S = 1/2 spins in a bipartite 
lattice one can in fact show a consequent reduction of the Hilbert space dimension from oo [as in 
linear spin-wave (LSW) theory] to ^^y- (as in MSW) which restores, up to logarithmic accuracy, 



the physical value of 2 

Takahashi's constraint imposes that (n)/(2S) < 1, in agreement with the kinematic constraint 
on the physical Hilbert space (even without explicit account of the projection operators on that 
space), and it guarantees the correctness of the truncations of high powers of n/(2S) that we intro- 
duced above. Finally, if the Hamiltonian is Z2 x U(l) symmetric because of a uniaxial anisotropy, 



as in the case of interest in this paper, one expects (s£) = 0. The constraint Eq. (10) elegantly 



restores this reflection symmetry of the ground state with respect to the quantization axis. 

B. Derivation of the self-consistent equations 

The correct spin wave description is found by minimizing the free energy T ' = E — TS , where 

S = '^[(n k + l)ln(n k + I) - n k lnn k ] (11) 
k 

is the entropy of a set of harmonic oscillators. Minimizing with respect to 9 k and u k under the 



constraint of Eq. (10) yields a set of self-consistent equations, 



with 



A, 

tanh2# fc = -^ (12) 

-Dfc 



A k = -^^Uj cos (Q ■r ij )G ij e ik - r ^ , (13a) 



N 
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where \x is the Lagrange multiplier for Eq. ( 10 ) corresponding to the chemical potential for changing 
the total magnetization. 



In Eqs. (13) we have abbreviated F^ = F(nj), and Gij = G(rij). Note that in the classical 
limit S — > oo one gets Gij,Fij S and Eqs. ( Jl~3| ) become analogous to their LSW counterparts. 
The spin-wave spectrum reads 



u k = JBl-Al. (14) 



Inserting Eq. (13) into Eq. (14) shows that a finite \x entails a gap at k = 0. This is to be seen in 
contrast to LSW theory where the spectrum always has a gapless Goldstone mode at k = 0. The 
correlators at the minimum take the form 

k 

Gi i = AT^Z^ cos(fc-r ij ) ^n fc + . (15b) 

At T = where = Vfc 7^ 0, one finds that [i vanishes. This implies also the disappearance of 
the gap at k = that may exist for finite temperature. A vanishing gap is a necessary requirement 
for the appearance of the Goldstone mode associated with magnetic LRO. It also enables Bose 
condensation of the DM bosons in the k = mode. Separating out the contribution of the zero 
mode, (a^, =0 afc = o) /N = (afc = oafc = o) /N = Mq (corresponding to the order parameter measuring 
the total spiraling magnetization in the quantization axis directions given by the ordering vector 
Q), one arrives at the zero-temperature equations 

Fij = M o + 7^7 £ — cos ( fc • r «) ' ( 16a ) 

Gij = M + -L V cos (fc • r -) , (16b) 
2iV Wfc 



and the constraint Eq. (10) becomes 



fc^O 



As mentioned above, the occupation of the zero mode Mo corresponds to a Bose-Einstein conden- 
sate of the DM bosons in the minimum of the dispersion relation. This condensate is depleted by 
interactions of the DM bosons. The larger this depletion, the more DM bosons reside at momenta 
different from zero, thus decreasing magnetic LRO. 
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C. Optimization of the ordering vector 



It is not a priori clear that the classical ordering vector Q cl correctly describes the LRO in the 
quantum system. To account for the competition between states with LRO at different ordering 
vectors Q we extend the MSW procedure by optimizing the free energy with respect to the ordering 
vector Q. This procedure, first introduced in Ref. |16j . significantly improves the predictions of 
MSW theory. It amounts to finding the best ordered reference state with in-plane ordering vector 
Q (spiral state) whose free energy is minimized not at the classical level, but including the effect 
of quantum fluctuations self-consistently within MSW theory. 

The minimization of T with respect to Q x and Q y yields two additional equations which must 
be added to the set of self-consistent equations, 





dQ. 



o, 



(18a) 



(id) 



d 



(18b) 



The values of Fij and Gij can now be calculated by solving self-consistently Eq. ( 18 ) together 



with Eqs. (10 12-15). At zero temperature, Eq. (10) and Eq. (15) have to be replaced by Eq. (17) 
and Eq. (16), respectively. Through Wick's theorem the knowledge of the quantities F^ and Gij 



allows for the computation of the expectation value of any observable. 



D. Spin stiffness 

The optimization of the ordering vector allows for a straightforward calculation of the spin stiff- 
ness. This additional information, complementary to the order parameter, helps us in identifying 
candidate regions for spin-liquid behavior. 

The MSW Ansatz always returns only a single one of all the possible ordering vectors Q° as the 
optimal ordering vector. However, if the true ground state is only short-range ordered, we might 
expect the Q-minimum to be relatively shallow, and that a slight change of the ordering vector 
barely affects the free energy T . This means that the order is not very stable against twists of the 
spin configuration. We quantify the curvature of the minimum by the spin stiffness tensor 

1 d 2 ^ 



Pa/3 



N dQadQp 



(19) 
Q=Q° 
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evaluated at the optimized ordering vector Q°. In particular we will determine the parallel spin 
stiffness pn = ^ Tr p = ^ (fa + p yy ) and the Gaussian spin stiffness 

T = det p . (20) 

The spin stiffness gives a measure of how stiff magnetic LRO order is with respect to distortions 
of the ordering vector and it provides a fundamental self-consistency check of our approach. In fact, 
finding a small spin stiffness casts doubt on the reliability of the spin-wave approach in describing 
such a strongly fluctuating state, and hence suggests that the true ground state might be quantum 
disordered. 

Since a change in Q affects the correlators and Gij, we must compute T self-consistently. 
After finding the optimal Q° by the self-consistent procedure described in the previous sections, 
we calculate jjT (Q x ,Qy) self-consistently for several ordering vectors Q = Q° + AQ and fit a 
quadratic form to the results. Since the minimum in the free energy can be very shallow, this 
procedure can be somewhat affected by numerical noise. As an approximation to the true spin 
stiffness, the partial spin stiffness /?^ rtial can be computed via the partial derivatives, i.e. without 
recalculating the self-consistent equations. It reads 

plT x - llJ^T^ (21) 



We define T partial analogously to T [Eq. (20)] as the determinant of the partial spin-stiffness tensor. 



The system can lower its energy by adjusting Fij and Gij to the new ordering vector, and therefore 
T is always smaller than T part . We will see later that in some cases the partial Gaussian spin 
stiffness T partial gives a good estimate of the real Gaussian spin stiffness T, but there are cases 
where it is considerably larger. 

In the following we present the zero-temperature phase diagram of the anisotropic triangular 
lattice. This paradigmatic example will show that the spin stiffness is an important quantity that 
provides a deeper insight into the order properties of the system. 



E. From spins to bosons 

As we mentioned in the introduction, in the 5 = 1/2 case the spin Hamiltonian is equivalent 
to that of infinitely repulsive bosons at half filling with frustrated hoppings. Hence it is important 
to match spin observables with their bosonic counterparts. Following the Dyson-Maleev or the 
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Holstein-Primakoff transformation, spin operator and physical (hardcore) b bosons obey the rela- 
tionship: = b\, = bi, Sf = Sf = b\bi — 1/2 [15] . Here the b operators obey anticommutation 
rules on site, {bi,b\} = 1. A non-zero magnetic order parameter Mo implies the appearance of off- 
diagonal LRO in the bosonic one-body density matrix, (b\bj) - — > Mq cosQ • ry. The ordering 
vector Q corresponds to the finite momentum at which condensation occurs. The condensed state 
in the spiral phase is characterized by a pattern of persistent currents forming a crystal of vortices, 
whose geometrically correlated structure is captured by the spin chirality (see below). Finally the 
parallel spin stiffness corresponds to the superfluid density of the bosons, p s = p\\/S. 

III. GROUND STATE PHASE DIAGRAM OF THE ANISOTROPIC TRIANGULAR 

LATTICE 

In this section we compute the ground-state phase diagram of the spatially anisotropic triangular 
lattice (SATL) with nearest-neighbor (NN) XY interactions. We consider a wide range of a = t2/t\, 
where t± denotes the bond strengths along the chains and t% the bond strengths along the diagonals 
(black and red bonds, respectively, in Fig. [I]). The parameter a interpolates between an ensemble 
of decoupled one-dimensional chains at a = 0, the isotropic triangular lattice at a = 1, and the 
square lattice for a — > oo. 

If we assume the spins to behave classically, the 2D-Neel order, present for a > 2, starts to 
continuously deform into spiral order at a < 2 [compare Fig. [2] (a)]. The spiral phase extends 
down to a = where the chains decouple. In a previous publication we presented the quantum 
mechanical phase diagram as predicted by projected entangled-pair states (PEPS) calculations 
[llj ; it is reproduced for convenience in Fig. [2](b). According to this, both the square lattice limit 
(a — > oo) and the most frustrated case, the isotropic triangular lattice (a = 1), display magnetic 
LRO. In the limit of uncoupled chains (a = 0) the system displays quasi-LRO with algebraically 
decaying correlations. However, similarly to what has been found in the Heisenberg model [24 , 
the system seems to feature spin-liquid phases with exponentially decaying correlations between 
different types of order or quasi-order. In Appendix A we provide a further spectral feature, coming 
from the exact diagonalization of a small cluster, which is consistent with the observation coming 
from PEPS calculations. Further distinct features of the quantum model are that the transition 
between 2D-Neel and spiral order is shifted by quantum fluctuations to considerably smaller values 
of a, and that the quasi-ordered ID-like state extends over a whole region of finite a in the phase 
diagram. 
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~ 0.4 ~ 0.6 ~ 1.2 ~ 1.4 t 2 /h 



Figure 2: (a) Classical ground-state phase diagram of the anisotropic triangular lattice with sketches of the 
ID state at a = 0, the spiral state at a = 1 and the 2D-Neel state for a > 2. (b) Quantum mechanical 
phase diagram from Ref. [11]. SL is short for spin liquid. 

In the following we compare predictions of MSW theory with these PEPS results and exact 
diagonalizations (ED). This will allow us to validate the reliability of the MSW method. The 
following system geometries are considered for the three different methods: 

• PEPS: a rhombic lattice of 20 x 20 = 400 spins with open boundary conditions and bond 
dimension D = 2. PEPS is a powerful numerical tool which goes beyond mean- field theory, 
but which, for small bond dimension D, only partially accounts for the entanglement in 
the ground state. This limitation becomes particularly serious close to quantum phase 
transitions. However, in Ref. [11] it was demonstrated that D = 2 is already accurate 
enough to effectively capture the physics of the system. 

• ED: Lanczos diagonalization of clusters of 24 and 30 spins (the latter is shown in Fig. [3]), 
again with open boundary conditions (necessary in order to allow for the accommodation of 
arbitrary ordering vectors). 

• MSW: rhombic lattices of 32 x 32 = 1024 and 64 x 64 = 4096 spins and in the infinite-lattice 
(thermodynamic) limit, under periodic boundary conditions. We find that at these lattice 
sizes all quantities have essentially reached the infinite lattice limit in most of parameter 
space. As it can be expected, the deviations from the thermodynamic limit are sizable only 
in the one- dimensional limit and at critical points. The thermodynamic limit is computed 
by replacing finite sums over the first Brillouin zone with integrals. 

In the triangular lattice with nearest-neighbor interactions one finds in the MSW formalism 
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Figure 3: Cluster of 30 spins for which we carried out exact diagonalizations (ED). The 24-spin system is 
equivalent, only with the top and bottom rows removed. The clusters are chosen for largest symmetry with 
respect to reflection on the axes and for a ratio of ti- (red) to ii-bonds (black) closest to 2. 



that Eq. (18b) gives Q y = and from Eq. (18a) we obtain the formula 

T n F 2 + f7 2 1 

(22) 



Q x = 2 arccos 



2FI +G2 



Here t\ = (1,0) and T2 = (l/2, v3/2) are the primitive lattice vectors (here and in the rest of the 
work we take the lattice spacing a equal to unity). For Fij = Gij = S, attained when S — > oo, this 
reduces to the ordering vector of the LSW theory (Q x l , Q~) = (2 arccos (—a/2) , 0) which coincides 
with the classical ordering vector. 

A. Breakdown regions for MSW theory 

As a first step in our analysis we investigate the parameter regions where LRO is to be expected, 
and the locations where MSW theory ceases to be applicable. To this end we first investigate if there 
appear imaginary modes in the dispersion relation, which would indicate instabilities. Afterwards 
we study the order parameter Mq and the spin stiffness. 

1. Imaginary frequencies and breakdown of convergence 

Convergence in the self-consistent equations of MSW theory with ordering vector optimization, 
Eqs. ( [T2}|T4| [l6| [Tt| [22] ), cannot be achieved in selected regions of the ground-state phase diagram, 
namely for a < 0.18 and for 1.35 < a < 1.66, as summarized in Fig. |4j (Interestingly, convergence is 
restored in the pure ID limit, a = 0, for which the theory formulates surprisingly good predictions.) 
This breakdown of convergence corresponds to the appearance of an imaginary part in the spin- 



wave frequencies, Eq. (14), signaling an instability of the ordered ground state. The breakdown of 
a self-consistent description of the system in terms of an ordered ground state is strongly suggestive 
of the presence of a quantum- disordered ground state in the exact behavior of the system. Hence, 
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one can interpret these parameter regions as candidates for the spin-liquid phases predicted by 
PEPS calculations [H] [compare Fig. ]2](b)]. 

■ , ■ J: 

0.18 1 1.35 1.66 2 a 

Figure 4: Regions in the phase diagram in which imaginary modes appear and convergence of the MSW 
equations breaks down (red). 



2. Order parameter and spin stiffness 

A fundamental indication for the validity of spin-wave theories is generally given by the order 
parameter Mq (Fig. [5]) and the spin stiffness (Fig. [6]). The influence of quantum fluctuations is 
strong where they are small, and the primary assumption that the system can be described by a 
semi-classical spin-wave state begins to falter in such a case. Since the MSW formalism only takes 
quantum fluctuations partially into account, a small order parameter and/ or spin stiffness also 
suggests that the true quantum ground state could be completely disordered. 

Interestingly in the regions of largest spatial isotropy of the interactions, i.e. around a = 1 
(isotropic triangular lattice) and at large a (isotropic square lattice), the order parameter Mq 
coincides with that of LSW theory. As we will see later, at the points a = 0, 1, oo the ordering vector 
found with the present MSW approach also exactly matches the classical one, due to symmetry. 

In the square lattice limit a — > oo, the order parameter attains the value Mq = 0.435 in the 
thermodynamic limit, which is very close to M = 0.437 as extrapolated from quantum Monte Carlo 
calculations |25j . For the spin stiffness Ref. [25] obtained p\\/a = 0.270; the MSW method returns 
the only slightly larger value p\\/a = 0.272. It appears that in this special case the main quantum 
corrections are correctly captured by MSW theory. The large values of the order parameter (around 
87% of the theoretical maximum) and of the spin stiffness support the assumption that the classical 
picture remains essentially valid in the large-a limit. 

The loss of LRO in the ID-limit (a — > 0) is reflected in the breakdown of the order parameter Mq, 
which occurs at a finite value of the inter-chain coupling, a ~ 0.18 (note that within LSW theory 
the order parameter vanishes only for a — > 0). This coincides with the appearance of imaginary 
spin-wave energies as discussed in the previous section. At small but finite a the spin stiffness p yy 
essentially vanishes, which is characteristic of a ID-like state that consists of effectively decoupled 
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Figure 5: Comparison of the MSW order parameter M [as defined before Eq. (16l] for different system 



sizes (numbers behind the labels). Also included is y/ S {Q) JN [Eq. ( 23 1] for ED (Lanczos) and PEPS 
computations, and the staggered magnetization M of LSW theory. A large value indicates strong LRO, 
with the theoretical maximum being 0.5. 



chains. This suggests that the physics becomes basically independent of the y-component of Q for 
a < 0.35. 

A single XY chain can be solved exactly by Bethe-Ansatz equations, and by use of twisted 
boundary conditions one can obtain the exact solution for the spin stiffness p xx = 1/vr ~ 0.318 
[26]. Our MSW result of 

Pxx ~ 0.308 lies surprisingly close. For one-dimensional models it is 
known that a non-zero spin stiffness is accompanied by quasi long-range correlations with power- 
law decay. The critical nature of the state in the ID-like phase reflects itself also in the fact that 
finite-size effects play an important role. 

The MSW order parameter can be compared with results from exact diagonalization (ED) and 
PEPS calculations. In both cases the Fourier transform of the spin-spin correlations, the static 
structure factor 



(23) 



allows to extract an ordering vector Q which maximizes S (k), as well as the order parameter M 
which is defined as M = \J S (Q) /N in the thermodynamic limit. In Fig. [H] we compare \J S (Q) /N 
for ED of systems of 30 and 24 spins, and with PEPS calculations on a 20 x 20 lattice. 

The comparison of ED and MSW results shows that MSW is quantitatively reliable in the Neel 
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Figure 6: (a) Gaussian spin stiffness T; (b) components of the spin stiffness tensor. The mixed second 
derivative p xy vanishes for symmetry reasons. The numbers behind the labels of the graphs give the consid- 
ered system size. The inset in (a) is a close-up of the region of small a. The curves labeled 'partial' were 



obtained by Eq. (21 1 



phase for a > 1.66. For smaller a values the comparison is more problematic: in particular, while 
ED and PEPS confirm the existence of an ordered spiral region for a around 1, the magnitude of 
the order parameter appears to be largely overestimated by MSW theory, which is not surprising 
considering the partial account of quantum fluctuations by the MSW approach. In particular, MSW 
theory produces the counterintuitive prediction that the frustrated spiral phase for 0.18 < a < 1.35 
has an order parameter which can be larger (around the isotropic a = 1 point) than that of the 
unfrustrated case of the square lattice (recovered for a — > oo). We observe that MSW theory is 
only moderately improving upon linear spin wave theory around the a = 1 point for what concerns 
the quantum fluctuations of the order parameter - in particular, its prediction for Mq essentially 
coincides with that of LSW for a = 1. 

In summary, from the MSW order parameter Mq we can derive a loss of LRO at a < 0.18, 
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and the spin stiffness suggests a strong weakening of inter-chain correlations already at a < 0.35. 
The spin-stiffness also decreases strongly upon approaching the parameter region 1.35 < a < 1.66. 
Together with the appearance of imaginary spin- wave frequencies for a < 0.18 and 1.35 < a < 1.66 
this strongly indicates the appearance of disordered phases in these regions. Due to its semi classical 
character the MSW Ansatz is not adapted to properly describe these phases, and we must resort to 
methods such as PEPS which take quantum fluctuations into account more completely. However, 
in the rest of parameter space magnetic LRO order seems to survive quantum fluctuations. 
In the next section we investigate in detail the nature of the ordered phases. 



B. Ordering vector, spin spin and chiral correlations 

In this section we introduce several observables which reveal the type of order appearing in the 



system, and which will be used in the discussions of sections III C and III D 

The ordering vector is a direct outcome of MSW theory, and can be extracted from the ED data 
by determining the position of the peak of the static structure factor, Eq. (23). Figure [7] displays 
the x-component of the ordering vector Q (with Q y = 0). Three limiting values are known. For 
a = intra-chain antiferromagnetic (Neel) order is described by Q = ttx. For a — > oo square-lattice 
Neel order is described by Q = 2irx. In the isotropic lattice (a = 1) the threefold symmetry forces 
the ordering vector to Q = (The ED and PEPS results deviate at a = 1 because the required 
threefold symmetry is broken by the shape of the simulation clusters, Fig. |5J) The importance 
of optimizing the ordering vector is apparent in Fig. [7] when comparing the MSW results to the 
classical (and LSW) curve. 

The spin-spin correlations of nearest neighbors shed further light on the order properties. We 
analyze them through the two-site total spin, 

K ij = ±{(S i + S j ) 2 ) = {S i -S j ) + l. (24) 

In Fig. [8] we plot it for nearest neighbors. This quantity vanishes if the spins are in a singlet, which 
is equivalent to perfect anticorrelation, takes the value | if they are uncorrelated, and the value 1 if 
the spins form a triplet, which means perfect correlation. For PEPS and ED we report the values 
of Kij averaged over the central spins, where boundary effects are minimal. 

Spiral phases carry not only a magnetic order parameter, but also a chiral order param- 
eter. In particular, a vector chirality [27 can be defined on an upwards pointing triangle 
with counter-clockwise labeled corners (i,j,k) as ka = [S{ x Sj + Sj x Sk + Sk X Si] z , 
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Figure 7: Comparison of the cc-component of the ordering vector, Q x , using ED (blue), PEPS (light green) 
and the MSW Ansatz (red). Also shown are the classical values (black). The numbers in the labels of 
the curves are the respective system sizes. The black circle marks the isotropic spiral ordering vector of 
Q x = 120° which occurs classically and within MSW theory at a = 1. 



and on a downwards pointing triangle with counter-clockwise labeled corners as Ky = 

[Si x Si + Si x Sj + Sj x Si] z . Long-range chirality correlations are defined as [28] 

tp- = ((«a - «v) ( K A' - «V')) ' ( 25 ) 

where the triangle pairs (A,V) and (A',V') share a T\ edge. In Fig. [9] we plot the average 
chirality correlation of the central plaquette with all other plaquettes, normalized to the theoretical 
maximum 4/9. The MSW data have been obtained by expanding the chiral correlation up to 
the fourth order in the boson operators, which is consistent with the truncation of the bosonic 
Hamiltonian Eq. (|5]) to the same order. Going to higher orders does not change the outcome in 
the regions where Mq is large, but can yield different results where Mq is small. In particular, the 
unphysical negative values attained by for small a are an artifact of this truncation. 

C. Transition from 2D-Neel order to spiral order 

An inspection of Figs. [7J [8] and [9] shows that the MSW formalism indeed reproduces the main 
features of the phase diagram of Fig. [2] (b) quite accurately. 

First of all, coming from the large-a limit, we observe that all methods (MSW theory, ED and 
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Figure 8: Nearest-neighbor correlation Kq^, where Xj — t\ (solid lines) and Xj = X2 (dashed lines), 
respectively, comparing ED (blue), PEPS (green), and MSW theory (orange). The black triangles are the 
MSW data for a one-dimensional chain of length N = 10 6 and the stars in magenta are the exact results for 
a linear chain in the thermodynamic limit derived by use of the Jordan-Wigner transform. 




Figure 9: Averaged chiral correlation normalized to the theoretical maximum 4/9 for ED (blue), PEPS 
(light green), and MSW theory (orange, dark green, and red). The black dotted line is the classical result 
and the black solid line is the classical chiral correlation that is obtained if for a given a the Q of the MSW 
calculation rather than Q cl is used. 
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PEPS) show a jump in the wavevector associated with the dominant correlations in the system, from 
Q = 2ttx, characteristic of the Neel phase, to a continuously varying Q, characteristic of a spiral 
phase. PEPS indicates a jump from dominant Neel correlations to dominant spiral correlations 
at a ~ 1.4. ED for 30 spins, shows a first order phase transition (with a sharp level crossing) 
between Neel and spiral correlations at the slightly larger value of a ~ 1.44. In the case of MSW 
theory, the jump in the ordering vector is realized when going across the breakdown region, namely 
when passing from a ~ 1.66 (which is the lower bound to the Neel phase within MSW theory) to 
a rj 1.35 (which represents the upper bound of the spiral phase). In particular, all three different 
approaches point to the fact that Neel order persists to much lower a than the classical value a = 2. 

The persistence of Neel order over a larger parameter region than in the classical case is remi- 
niscent of what is observed in other models. Indeed, quantum fluctuations generally stabilize states 
where spins are ordered collinearly (see e.g. Refs. [29] [30]), a property that is reproduced by the 
MSW Ansatz with ordering vector optimization. The mechanism behind it is strongly connected 
to order- by-disorder phenomena [30 . 

The abrupt transition from a phase with dominant Neel correlations to a phase with dominant 
spiral correlations is confirmed by the spin-spin correlations as displayed in Fig. [8j Anti-correlation 
along the strong T2-bonds and correlation along the weak Ti-bonds are characteristic of a 2D-Neel 
ordered state; these correlations decrease rapidly for a < a cnt outside of the Neel-ordered phase. 
The change of the type of order is further supported by the overlap KV'alV'oo)! of the new ground 



state with the 2D-Neel ordered state of a = oo, which we plot in Fig. 10 for the ED: above the 
phase transition it still attains a finite and quite large value, while it vanishes identically below 
the phase transition. Finally, the onset of strong chiral correlations shows that the new phase is 
indeed a spirally correlated one (Fig. [9J . 

The breakdown region of MSW theory, 1.35 < a < 1.66, is strongly suggestive of the loss of 
magnetic LRO, corresponding to a spin-liquid state. This region of disordered behavior is only 
roughly consistent with that indicated by PEPS calculations [11] for the appearance of a short- 
ranged spin-liquid phase, namely 1.2 < a < 1.4. Nonetheless it is tempting to associate the 
breakdown of MSW theory to this quantum-disordered phase. 

D. Persistence of ID quasi-LRO up to finite inter-chain couplings 



In MSW theory, at a ~ 0.18 the order parameter Mq breaks down. This is an indication of the 
transition to a phase without magnetic LRO such as the one reproduced in the one-dimensional limit 
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Figure 10: Overlap of the ground state at a with the 2D-Neel ground state of a = oo, |^a;=oo)j an d with 
the six-dimensional subspace that corresponds to the ground state of a = 0, respectively (ED, 30 spins). 

a — > 0. For a > 0.18 the weak a-dependence of the ordering vector, the spin-spin correlations and 
the energy strongly indicate the persistence of a phase with ID quasi-LRO properties to finite inter- 
chain couplings. Furthermore, the spin stiffness component p yy practically vanishes for a < 0.35 
which suggests that, if we took quantum fluctuations into account in a more accurate way than in 
MSW theory, spiral (2D) order would be probably lost below a < 0.35. 

In the limit a — * we find that the intra-chain spin-spin correlation Ki^ +Tl = K Tl ~ 0.355 
lies close to the exact result that can be obtained by use of the Jordan- Wigner transform, K Tl = 
— - — ^2 + f — 0.330. The inter-chain spin-spin correlations on the other hand vanish, which is 
equivalent to Ki^ +T2 = K T2 = |. 

We find that the ordering vector of the MSW theory compares very well to the one computed 
by ED in the entire range of a (Fig. [7|. Especially the very weak dependence of Q on a near the 
ID-limit is found in both approaches, contrary to classical and LSW theories which exhibit linear 
dependences on a. This means again that quantum fluctuations stabilize collinear order within 



the chains. This is confirmed by ED: the overlap y^ i=1 |(V ; a|V'o)| 01 t ne ground state with the 
subspace spanned by the six-fold degenerate [3] ground-states of a = remains very large (almost 



80 percent) up to a ~ 0.5 (Fig. 10). Moreover, the small chiral correlations in both PEPS and ED 



show that for small a there is no spiral long-range order. 



23 



E. Momentum distribution of the hardcore bosons 



After characterizing the zero-temperature phase diagram of the spin model in the previous sec- 
tions, we now wish to make contact with the cold-atom implementation of such a model via hardcore 
bosons. The most common observable in cold- atom experiments is the momentum distribution [I], 



which exactly corresponds to the static structure factor for S = 1/2 spins (23) 



,l >> (*) = Jf E e-*^-"> (b\b,) = S (*) (26) 



via the spin-to-hardcore-boson mapping described in Sec. II E Figure 11 shows the MSW prediction 
for the momentum distribution at various a values, spanning all the condensation regimes of the 
bosons at zero temperature. At a = (not shown) the system displays quasi-condensation at finite 
momenta along the uncoupled chains, resulting in vertical ridges at Q x = ±7r in the momentum 
distribution. These ridges corrugate as the interchain coupling increases, and true condensation 
peaks emerge in reciprocal space, corresponding to a condensate state which supports a crystalline 
vorticity pattern. For a = 1 these peaks are located at the six corners of the first Brillouin zone. 
For a < 1 the peaks are elongated in the y direction, while for a > 1 they are elongated in 
the x direction, witnessing the spatial anisotropy of the lattice. This situation persists up to the 
breakdown of MSW theory at a = 1.35; after recovery of the theory at a = 1.66, the momentum 
distribution shows condensation at the four corners of the Brillouin zone of a (deformed) square 
lattice, defined by the dominant diagonal bonds of the spatially anisotropic triangular lattice. 

The peak height (normalized to the number of sites) is given by the square of the order parameter 
M . 



F. Discussion 



Here we summarize the main features of the zero-temperature phase diagram obtained via MSW 
theory with Q-vector optimization. The region where the system behaves like a set of essentially 
decoupled chains with ID quasi-order is extended to considerable inter-chain interactions. The 
order parameter indicates that the point where inter-chain correlations set in occurs at a ~ 0.18; 
the spin stiffness suggests that an effective decoupling of the chains may even persist up to a ~ 0.35. 
At larger a the system crosses over to a spirally ordered phase that persists up to a ~ 1.35, where 
MSW breaks down, suggesting a quantum disordered ground state. At a ~ 1.66 MSW theory finds 
again a self-consistent solution, this time corresponding to a 2D-Neel state. 
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Figure 11: Momentum distribution of a half-filled gas of frustrated hardcore bosons on the spatially 
anisotropic triangular lattice [corresponding to the spin static structure factor S(k)], for various values 
of the spatial anisotropy a. The data (on a logarithmic color scale, scaled to the number of sites) result 
from a MSW calculation on an 18 x 18 lattice. The black hexagon marks the first Brillouin zone and the 
black cross its origin. 



These results are mostly consistent with the PEPS phase diagram of Fig. @ (b). Especially 
the persistence of ID behavior to surprisingly large values of a, the fact that the long-range 
ordered spiral phase survives quantum fluctuations, and the extension of 2D-Neel LRO down to 
much smaller values of a than in the classical equivalent are reproduced. However, there are some 
deviations, which are generally to be expected from a spin- wave approach. The range of the ordered 
phases appears to be somewhat overestimated by MSW theory. Furthermore the gapped spin-liquid 
phases are not faithfully described: while the breakdown of MSW theory for 1.35 < a < 1.66 
suggests a disordered ground state, the gapped spin liquid in the region 0.4 < a < 0.6, suggested 
by the PEPS data pT], is not observed in the MSW results. Still the proposed phase diagram is 
significantly improved with respect to the LSW phase diagram, which follows the classical phase 
diagram too closely. 

Two improvements have proven to be crucial: First, the minimization of the free energy with 
respect to Q in the self-consistent equations has enabled us to describe the considerable shift of 
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the ordering vector to a surprisingly satisfactory level. Second, the investigation of the Gaussian 
spin stiffness serves as a measure of the actual stiffness of the long-range ordered phase. As such 
it allows us to detect regions where spin-liquid behavior may appear in the true quantum ground 
state. 

IV. FINITE TEMPERATURE PHASE DIAGRAM OF THE ANISOTROPIC 

TRIANGULAR LATTICE 

We now investigate how the phase diagram of the antiferromagnetic nearest-neighbor (NN) XY 
Hamiltonian on the anisotropic triangular lattice extends to finite temperatures, making use of 
the modified spin-wave (MSW) theory with ordering-vector optimization. All calculations in this 
section are carried out in the thermodynamic limit. 

At finite temperatures, continuous symmetries cannot be spontaneously broken in two dimen- 
sions |31[ 132] . In the XY model, instead of long-range order (LRO), one finds quasi-LRO at finite 
but low temperature. At the Berezinskii-Kosterlitz-Thouless (BKT) temperature Tbkt the system 
undergoes a topological phase transition from quasi-LRO to an exponential decay of correlations, 
involving the unbinding of vortex-antivortex pairs |33[ IM) [35] . The existence of a BKT transition 
in the XY model must be seen in contrast to the Heisenberg model, where Kosterlitz and Thou- 
less showed that vortex excitations are not topologically stable [35 , a fact which precludes the 
possibility of a BKT transition. 

The possibility of observing the BKT transition is a particular advantage of MSW theory. 
The BKT phase with algebraic order is generally predicted by linear spin-wave (LSW) theory to 
remain stable at arbitrary temperatures. The non-linearities contained in MSW theory allow the 
disruption of quasi-LRO and the transition to the short-range-ordered (SRO) phase. However, 
vortex-antivortex excitations are not explicitly present in the theory, and therefore in principle 
Tbkt cannot be accurately estimated. 

A. Spin spin correlations 

An important observable for the analysis of a temperature-dependent phase diagram is the 
two-point correlation function 

c ij = (SfS* + SfS]) I cos (Q • nj ) = \ (i* + G%) . (27) 
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Figure 12: Correlation function C mT , along the chains in the thermodynamic limit, (a) for a = 0.7, and (b) 
for a = 100 for several temperatures. The normalized temperatures T/(t\ + It?) considered (for lines from 
top to bottom) are (a) T/(t 1 + 2t 2 ) = 0.064 . . . 0.078 in steps of 0.002, and (b) T/(tt + 2t 2 ) = 0.068 . . . 0.164 
in steps of 0.016. The line closest to Tbkt calculated from MSW data is shown bold. Note the difference in 
the maximal m between (a) and (b). 



In our analysis we focus on C mTl and C mT2 where m is a positive integer, and t\ = (1,0) 
and T2 = (l/2, y/E/2\ are the lattice vectors. The behavior of C mTl captures the intra-chain 
correlations, while that of C mT2 describes inter-chain correlations. 

In order to locate the BKT transition we calculate the residual sum of squares R = 
Ylm [Cmru — f( mT i,2)] 2 for two trial functions, an exponential f(r) = Ae~ r ^, where £ is the 
correlation length, and an algebraic fit f(r) = A/r 11 . We fit these functions to the correlations of 
the central spin with sites which are m = 3 ... 15 lattice spacings apart. The lower limit to the 
fit region is necessary because the trial functions are only valid for the long-distance part of the 
correlations, while the upper limit is chosen by us because computation time for the correlations 
increases considerably with the distance between the spins. We identify a BKT transition with 
a point where the residual sum of squares R of the exponential fit becomes equal to that of the 
algebraic fit. This has to be understood as only a rough estimate of the transition temperature. 
When giving explicit values of transition temperatures we take the average of the values obtained 
from fits to C mXl and C mX2 . 

Figure [l2] shows representative log-log plots of the correlation function C mTl at a = tijt\ = 0.7 
and a = 100 (where the ground states show spiral and 2D-Neel order, respectively) for several 
temperatures. In these plots algebraically decaying correlations correspond to straight lines. For 
a = 0.7 a clear transition from algebraic to exponential decay at the computed BKT temperature 
can be observed, a behavior which is found in the entire parameter range of the spiral phase. On 
the contrary, at a = 100 we cannot find such a clear transition. Rather, the curves acquire a 
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curvature in a fairly continuous way, which makes it difficult to pinpoint the transition. In order 
to check how the suggested BKT line changes when taking correlations to more distant spins into 
account, we computed spin-spin correlations for a. = 100 up to distances of 64 lattice spacings. 
This yields a transition temperature of Tbkt/^i +2i2) = 0.134, which is approximately 15% lower 
than what is obtained if distances of only up to 15 lattice spacings are considered. In light of the 
approximate nature of MSW theory, we find that this level of precision is satisfactory. 

Another useful observable is represented by the gap A = Afc = o of the spin- wave dispersion, 
which is intrinsically connected to the two-point correlations. The gap is directly imposed by the 



chemical potential \x [see Eqs. (13) and (14)], and its magnitude determines the rapidity of the 
decay of correlations. A finite gap leads to exponentially decaying correlations while a vanishing 
one entails power- law correlations. Hence in principle the onset of a gap at finite temperatures 
corresponds to the occurrence of a BKT transition. 

In reality, the thermal onset of a gap we observed within MSW is typically very gradual, and a 
clear identification of the transition point via the gap is generally problematic. This observation 
can be understood on the basis of a well-known fact: the chemical potential of the half-filled DM 
boson gas, which determines the existence of a gap, cannot vanish at finite temperature because 
of the absence of Bose-Einstein condensation in two dimensions. As a consequence, we find a 
finite gap at any finite temperature, which means that the correlations decay exponentially at 
long distances. This suggests that, strictly speaking, MSW theory is not able to describe the BKT 
transition. However, for temperatures much lower than the BKT transition (estimated as explained 
above) the gap is very small, being below our numerical precision. For all practical purposes such 
a small gap entails a decay of correlations which is not distinguishible from an algebraical decay. 
Moreover in a selected region of the phase diagram (corresponding to the spiral phase) the gap is 
seen to increase drastically around the estimated BKT transition temperature, and correspondingly 
the correlation function is seen to decay much more rapidly above that temperature. Hence we 
conclude that MSW theory still accounts for one of the most salient features of the BKT transition, 
namely a discontinuous behavior of correlations as the temperature is increased. 

Finally we can extract from the correlations the temperature at which the MSW formalism 
breaks down. It is characterized by the complete loss of all correlations, even to the nearest 
neighbor. This behavior, occurring at temperatures of the order of the coupling strength, is clearly 
an artifact of the method, since in real systems the complete loss of correlations occurs only at 
extremely large temperatures where spin-spin interactions become negligible. 
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Figure 13: Schematical temperature-dependent phase diagram of the XY SATL. The different regions are 
listed along with their main characteristics in table [T] Horizontal lines mark lD-Nel order; diagonal lines 
spiral order; and cross-hatches 2D-Nel order. 



Phase 


Q x 


decay of correlations 


(A) ID-like SRO 


7T 


intra-chain: exponential 






inter-chain: uncorrelated 


(B) Spiral quasi-LRO 


7r < Q x < 2ir 


algebraic 


(C) Spiral SRO 


7r < Q x < 2n 


exponential 


(D) 2D-Neel quasi-LRO 


Q x = 2tt 


algebraic 


(E) 2D-Neel SRO 


Q x = 2tt 


exponential 


(F) Unstable (imaginary modes) 

(G) Breakdown of theory 




no correlations 



Table I: Parameter regions found in the finite-temperature phase diagram of Fig. |13| We distinguish mainly 
between phases with quasi-long-range order (quasi-LRO), i.e. algebraic decay of correlations, and phases 
with short-range order (SRO), i.e. exponential decay of correlations. Moreover, two regions are listed where 
the MSW formalism ceases to be applicable [(F) and (G)]. 

B. The phase diagram 



In this section we present the finite-temperature phase diagram of the anisotropic triangular 
lattice model obtained via the MSW method with ordering vector optimization [JS] . We derive it 
from the observables introduced in the previous section. For reference we first present a summa- 



rizing sketch of the phase diagram in Fig. 13 which introduces the phase labels we will refer to in 
the following discussion. Table [I] lists the main properties of these phases. 

First, at small a, there is a phase with properties similar to the algebraic lD-Neel-like state 
found at T = but with exponential decay of intra-chain correlations (phase A). Further, we 
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generally find that the phase diagram contains two quasi-LRO regions: a region at intermediate a 
corresponding to spiral quasi-LRO (phase B), and another region at large a which is characterized 
by Neel quasi-LRO (phase D). These phases undergo BKT transitions to similar phases with 
short-range order (SRO), phases C and E, respectively. Moreover, between them lies a region 
where imaginary frequencies occur in the spin-wave dispersion relation, which can be interpreted 
as an indication for an extremely short-range-ordered phase (phase F). This general structure of 
the phase diagram is supported by all the observables we investigate. At large T the MSW method 



breaks down (see sec. IV A) and therefore does not allow for any interpretation in that domain 
(region G). 

It is also important to note that our calculations cease to converge properly for too low tempera- 
tures, when the chemical potential becomes smaller than the accuracy of our numerical integrations. 
Depending on the region of the phase diagram the lowest temperatures for which appropriate re- 
sults could be derived vary from less than one-tenth of a percent to several percent of the coupling 



strengths. This pathology is not observed at T = (as calculated in section III) because an exact 



vanishing of the chemical potential allows the special treatment of the zero-mode as captured in 



Eqs. (16) and (17). Except for some points, we typically calculated down to T/(t± + 2*2) = 0.025. 
Since the bond strengths are the only energy scales in the problem, it seems a reasonable assump- 
tion that our finite temperature results can be analytically continued down to T = + without 
encountering discontinuities (except possibly at exactly T = 0, where — in contrast to any finite 
temperature — Bose-Einstein condensation of the DM bosons becomes possible). Nevertheless, 
this issue should be kept in mind in the following analysis. 



The breakdown of the calculations for too low temperature can be clearly seen in Fig. 14 which 
displays the phase diagrams obtained from several observables. 

A natural starting point for a thorough analysis of the temperature dependent phase diagram 
is given by the respective ground state phases. Proceeding from small to large a, we divide the 
analysis in sections corresponding to different ground state behavior. 

1. ID-like phase (phase A) 

The one-dimensional quasi-ordered ground state phase for which we found strong indications 
below a ~ 0.18 becomes a short-range-ordered phase at finite temperature (phase A). It is char- 
acterized by vanishing correlations between neighboring chains already at zero temperature. The 
finite gap leads to exponentially decaying intra-chain correlations for all T. This is consistent with 
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Figure 14: Linear color plots in dependence of a and Tf (ti + 2t%) of (a) the x-component of the ordering 
vector Q X1 (b) the intra-chain correlation K Tl , (c) the inter-chain correlation K T2 , (d) the partial Gaussian 
spin stiffness T partlal , and (e) and (f) the partial spin stiffnesses p%% rtl£d and p^^* , respectively. The 
mixed component of the spin stiffness p!^ rtlal vanishes for symmetry reasons. Figures (g) and (h) show the 
entropy and the gap A. The points mark the BKT transition for T\ — (1,0) (red), the BKT transition 
for T2 = (l/2, V3/2) (orange), the break-down temperature (blue), and the temperature where inter-chain 
correlations disappear (yellow), all computed through the two-point correlations C mxi 2 . In the grey region 
imaginary frequencies appear in the spin- wave dispersion. 
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the expected finite-temperature behavior above a ground state with quasi-LRO. 

The assumption that this low-a phase really describes decoupled chains is reinforced by the 
component p?^ rtial of the spin stiffness, which vanishes in this region, and by the ordering vector 
that takes on the ID value (ir,0), similar to the equivalent of the ground state phase diagram. 
Moreover, neighboring spins on different chains are uncorrelated whereas nearest-neighbors on the 
same chain are anti-correlated. It is remarkable that this phase is preferred over the quasi-ordered 



spiral phase with rising temperatures. In section HID we have seen that quantum fluctuations 
stabilize lD-Neel quasi-order. The same mechanism is at work here: collinear spin correlations are 
stabilized by fluctuations, in this case thermal ones. 

Note that in the ID-like phase A the inter- and intra-chain correlations behave completely 
differently. In the rest of the phase diagram they follow one and the same pattern, since in a truly 
two-dimensional structure the correlations in one direction typically cannot disappear without 
affecting the correlations in the other one. 



2. Spiral phases (phases B and C) 

At intermediate inter-chain couplings 0.18 < a < 1.35 and low temperature we find a spiral 
phase with magnetic quasi-LRO (phase B). It can be seen as the finite-temperature continuation of 
the spirally ordered ground state phase. At larger temperatures a BKT transition to a phase with 
a spiral ordering vector but with an exponential decay of correlations occurs (phase C). Within 
our resolution of the phase diagram it seems that phase C disappears on the large-a side of the 
B-phase dome, and that phase B is delimited at large a by the F region where imaginary spin-wave 
frequencies appear. Furthermore, on the low-a side phase C becomes extremely narrow and is 
almost immediately followed by a transition to phase A described in the previous section. The 
broadest extent in temperature of C is around a = 1. 

At the isotropic point a = 1 the BKT transition from B to C is approximately located at 
Tbkt/(^i + 2^2) = 0.0836. Quantum effects lower the transition temperature considerably from the 
classical value Ibkt/C^i + 2t2) = 0.165 found by classical Monte Carlo simulations [36]. A pure- 
quantum self-consistent harmonic approximation, developed in Ref. [37 , gives Ibkt/(^i + ^2) = 
0.0625. The fact that MSW theory gives a significantly higher estimate is not surprising given that 
Ref. [37] takes vortex-antivortex excitations explicitly into account while MSW theory does not. 

We also find that in the same domain of large frustration (a close to 1) where spiral quasi-order 
is most stable, the breakdown of the theory occurs at a lower critical temperature. 
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We note also a strong drop of the transition temperatures around a ~ 0.4. In Fig. 13 this is 
marked by a dashed line which separates phase B from a phase B' with similar properties. We 
believe that this behavior is a numerical artifact and that in fact B and B' are one and the same 
phase. 



3. Spin-liquid candidate region (phase F) 

At the high-a side of the spiral phases we find an extended region where the spin- wave dispersion 
acquires imaginary modes. This means that MSW theory predicts an instability here. The width of 
this region in a stays approximately constant but it moves to smaller a with increasing temperature, 
leaving space to the collinear short-range-ordered phase (E). The region F extrapolates well down 
to the suspected spin-liquid phase between 1.35 < a < 1.66 at T = 0. Given that MSW is seen to 
break down at a putative spin-liquid phase at T = due to its lack of order, a fortiori one can 
expect MSW to break down in the same parameter range at finite temperatures, because at finite 
T the theory would be required to describe not only the ground state but also the excitations on 
top of it. 

Note also that the spin-stiffness decreases upon approaching this region, which could be inter- 
preted as a precursor of a short-range-ordered phase. 



4-. 2D-Neel states (phases D and E) 

As expected from BKT theory, when going to finite temperatures the 2D-Neel ground state 
first changes into a low-T quasi-long-range ordered phase (phase D), which at a temperature Tbkt 
undergoes a transition into a high-T short-range-ordered phase (phase E). Both are characterized 
by an ordering vector at the 2D-Neel value Q = (2-7r,0). Furthermore neighboring spins which 
share a diagonal bond are strongly anticorrelated whereas neighboring spins which lie on the same 
chain are positively correlated. 

The square XY lattice, which is reached as a = £2/^1 — ¥ oo, has been extensively studied in 
the past. The classical BKT-temperature T^ KT /t2 = 0.695, which has been calculated by use 
of classical Monte Carlo simulations [38 a , is significantly lowered in the quantum limit to around 
Ibkt/^ ~ 0.35 (Quantum Monte Carlo calculations [391 SO]). Our MSW results yield a BKT 
temperature of Tbkt/^2 ~ 0.27 at a = 100, where the system has practically reached the square 
lattice limit [46] . Once again the disagreement with the MSW result is not surprising, given that 
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this theory does not account properly for vortex-antivortex excitations. In particular the BKT line 



for a > 1.6 is not very distinct, due to the location problems mentioned in Sec. IV A Therefore its 
quantitative value should be interpreted with caution, especially in this parameter range. However, 
the qualitative behavior of the phase diagram seems to be described correctly. 

In the following section we turn to observables which show more clearly how order in the different 
phases persists at finite temperatures. 

C. Observables distinguishing between LRO and SRO 

Here we focus on observables which help to distinguish between different types of regimes (i.e. 
quasi-long-range order or short-range order), namely the partial Gaussian spin stiffness T partial , 
the gap A, the entropy, and the occupation of the zero mode rik=o- 

1. Entropy, spin stiffness, and gap 



The entropy [Eq. (11)] shows behavior consistent with the quasi-ordered character of the low- 



temperature 2D-Neel phase D and the spiral phase B in the sense that it is smaller in phases with 



stronger order as can be seen in Fig. 14 (g). Correspondingly, in phases B and D T partial is large 
[Fig. 14 (d)] and A is very small [Fig. 14 (h)]. 

At the BKT transition no sharp change occurs in these observables as may be expected. How- 
ever, the contour lines of the spin stiffness, the entropy, and the gap, all seem to be consistent with 
the shape of the Tbkt curve. 



We report the gap A for two representative values of a in Figs. 15 (a) and (b). It evolves 
smoothly through the BKT transition in the Neel phase, whereas in the spiral phase it is seen to 
display a sharp increase right above the BKT transition up to the transition to phase A. On the 
contrary, for all phases a sharp increase of A (accompanied by a sharp drop of "f partial ) can be 
discerned at the breakdown temperature where correlations are completely lost. 

In phase A (better visible for smaller values a which are not shown) we find that the gap is 
almost a linear function of temperature up to very close to the breakdown temperature. This is 
typical of critical systems where the temperature introduces the only energy scale, which is reflected 
in the gap. 
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Figure 15: Gap A [(a) and (b)] and occupation of the zero mode iik=o [(c) and (d)] for two values of 
a. The vertical lines denote the critical temperatures of the BKT transition calculated via the decay of 
correlations in direction T\ = (1, 0) (red), the corresponding BKT transition for x 2 = (1/2, v / 3/2) (orange), 
the break-down temperature (blue), and the temperature where inter-chain correlations disappear (black) 
[in (b) and (d) the last two markers overlap], all computed through the two-point correlations C mTl2 . 
Capital letters refer to the phases of Fig. [13] 

D. Occupation of the zero-mode 



The occupation of the zero-mode rifc = o [Figs. 15 (c) and (d)] gives an insightful measure of the 



strength of correlations. Due to constraint (10) a large population of the zero mode entails smaller 



population of excited modes and therefore leads to stronger correlations. 

Similar to what was seen for the gap in the previous section, at the transition between phases 
C and A rafc=o drops from very large values to something of the order of 1. Afterwards it changes 
only slightly with T. It is useful to remember that the average mode occupation, Y2k n k/N , equals 



S by virtue of constraint (10). This means that in the ID-like short-range-ordered phase A 71^=0 



is relatively small but still larger than the occupation of the other modes. 

The behavior of rifc=o is different for the 2D-Neel phase. At the BKT line that we extracted 
from the analysis of the two-point correlation functions C m Tl and C m T2 , ny,=o decreases strongly 
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but smoothly. Up to the breakdown point its values are still several times larger than in the ID- 
like phase, however. This supports our identification of the BKT transition; but the smoothness 



of nfc = o also shows the reason why the observables of section IV C could not point out a sharp 
transition. 



E. Discussion 



The finite-temperature phase diagram is observed to be a natural extension of the ground- 
state phase diagram. We find that zero-temperature long-range ordered phases are reflected in 
finite-temperature phases with quasi-LRO, while phases with quasi-LRO at T = turn into short- 
range-ordered phases at any finite temperature. When temperatures are at or below a few percent 
of the coupling strengths, the main characteristics of the ground state phase diagram are retrieved, 
with a short-range ID-like phase (A), and two quasi-ordered phases [one with spiral properties 
near the isotropic triangular limit (B) and one with 2D-Neel-like characteristics at large values 
of a (D)] which are separated by a potential spin liquid (F). This last phase was identified by 
(i) the breakdown of MSW theory, which indicates that the assumption of an underlying ordered 
state is invalid, and (ii) the lowering of the spin stiffness as this phase is approached. In [A] we 
show further evidence for spin-liquid behavior in this phase at T=0 gleaned from very limited 
exact-diagonalization results. 

We have given a rough estimate for Tbkt over the entire range of anisotropies. We find agree- 
ment in the rough magnitude of Tbkt at points where estimations of the BKT temperatures 
computed by other methods exist. In our results the BKT transition is more clearly visible in the 
correlations for the spiral phase than for the 2D-Neel phase. 



V. CONCLUSIONS 



We have extended Takahashi's modified spin-wave theory by an optimization of the ordering 
vector, which allows to account for order that deviates from the classical one. 

We have used this method to calculate the ground state phase diagram of the spatially 
anisotropic triangular lattice with S = 1/2 spins and XY interactions. We found the expan- 
sion of a quasi-ordered ID-like phase to finite inter-chain couplings, a spiral phase, and a 2D-Neel 
phase. At the transition between the latter two the breakdown of MSW theory indicates the loss 
of LRO. 
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We have extended this phase diagram to finite temperatures and computed Berezinskii- 
Kosterlitz-Thouless transitions, although the results are to be interpreted only semi-quantitatively 
because MSW theory does not explicitly account for vortex-antivortex excitations. We find that 
the ground state phases clearly imprint their properties on the finite temperature phase diagram 
at low temperatures, with long-range ordered phases being replaced by quasi-ordered ones and 
quasi-ordered phases by short-range ordered ones. 

Qualitative and even quantitative agreement with PEPS and ED calculations was found in 
the regions where magnetic LRO is to be expected. In particular it has been shown that MSW 
theory with ordering vector optimization is able to account satisfactorily for the main quantum 
corrections to the ordering vector. Furthermore, our calculations show that the spin stiffness is a 
useful observable for finding candidate regions for spin-liquid behavior in the ground state. Indeed 
the breakdown of MSW theory, or the very weak stiffness of the magnetic order it predicts, can be 
used as indications of the absence of long-range order in the exact ground state of the model. 

We find two main recurrent features for strongly frustrated quantum-magnets in two dimensions: 
collinear order is considerably stabilized by quantum and/ or thermal fluctuations against spiral 
order, and ordered or quasi-ordered phases characterized by different forms of order (collinear vs. 
spiral) do not continuously connect to each other, but they rather seem to be separated by quantum 
disordered phases. While MSW theory cannot determine the properties of such disordered phases, 
it provides a fast and clear method for finding candidates of disordered phases. This method can 
therefore serve as a guide in our search for interesting quantum-mechanical lattice models which 
require an experimental quantum simulator for further study of their phase diagram. 
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Appendix A: SIGNATURES OF ORDERING AND SPIN-LIQUID BEHAVIOR IN THE 
EXACT DIAGONALIZATION SPECTRA OF A SMALL CLUSTER 

We present here exact diagonalization data for the spectrum of the S =1/2 XY antiferromagnet 
on the spatially anisotropic triangular lattice. We perform our calculations on a 24-spin cluster 
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Figure 16: Spectra of the 5=1/2 XY antiferromagnet on the spatially anisotropic triangular lattice from 
exact diagonalization on a 24-site cluster, for various values of the lattice anisotropy a. The lower-right 
panel shows the average maximal level spacing (see text) as a function of a. 

with the geometry depicted in Fig. [3} The system Hamiltonian, Eq. (j2|, commutes with the total 
magnetization along the z axis, S£ , so that excited states can be classified on the basis of this 
quantum number |47] . 



Figure 16 shows the excitation energies of the first excited states in each 5*°* sector (up to 
gtot _ -q) -^jth re spect to the minimum energy in each sector, Eo(Sl ot ) [£'o(£'z 0t = 0) corresponds 
to the ground state energy]. Upon varying a we observe a significant evolution of the low-energy 
spectrum of the system, which points at the widely different regimes explored by the system. 
In particular, in the spirally and Neel ordered phases - exemplified in Fig. [16] by a = 1 and 
a = 2, respectively - we observe that in each S 1 * * sector there are a few states lying close to 
the minimum energy one, and separated from the other excited states by a large gap. According 
to a standard 'tower-of-states' argument j5T], these low-lying states are expected to collapse to 
the ground state in the thermodynamic limit, giving rise to degenerate superpositions of all S** ot 
sectors, each breaking the U(l) rotational symmetry of the Hamiltonian and displaying spiral or 
Neel order. The higher-energy states will instead reproduce the true excitation spectrum in the 
thermodynamic limit. 

This tower-of-states feature is on the contrary absent in other regions of the phase diagram, 
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in which the energy levels in each 5* ot sector are more homogeneously spaced. The absence of a 
low-lying multiplet of states separated from the higher energy states by a large gap is observed in 
models whose ground state is generally considered to be a spin liquid [42] . We therefore introduce an 
observable aimed at quantifying the extent to which the spectrum exhibits the expected features in 
presence of spontaneous symmetry breaking in the thermodynamic limit. We consider the average 
maximal level spacing A max , defined as 

N s 



,„,, ^— £ max[£ m (Sr)-£i(^ ot )] , (Al) 



A """ N s + 



namely A max is the maximal level spacing in each S z sector, averaged over the Ns + 1 = 12 
sectors considered. The maximal level spacing is extracted by considering the lowest 10 levels 
Ei(Sl ot ), which captures the behavior of the low-energy part of the spectrum. The above quantity 
is chosen so that it will be maximal in presence of a large separation between the low-lying tower 
of states and the higher-energy spectrum, while it will be minimal in presence of homogeneously 
spaced levels in each sector. 



When plotting A max as a function of a, as shown in Fig. 16 we observe two very pronounced 
relative minima, at a ~ 0.6 and a ~ 1.4. Remarkably these two minima correspond to the two 
regions in parameter space where PEPS calculations predict the occurrence of a spin-liquid phase 
jllj (compare Fig. [2]). Hence the lack of the tower-of-states feature in the spectra of a small cluster 
is consistent with the PEPS prediction. 
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